Whiteline: Alpha

home *** CD-ROM | disk | FTP | other *** search

/ Whiteline: Alpha / Whiteline Alpha.iso / progtool / modula2 / complex / cmplx_16.txt < prev next >

Wrap

Text File | 1994-09-22 | 18.2 KB | 361 lines

(C) 1989 Thomas Groß, Lynarstraße 13, D 1 Berlin 65 ____________________________________________________________ Anleitung zu CMPLX_16.FOR __________________________________ CMPLX_16.FOR ist eine Programmsammlung für doppelt-genaue komplexe Ausführung der INTRINSIC-Funktionen von Fortran. Darüber hinaus enthält sie einige in Fortran 77 nicht vorgesehene Funktionen. Es gibt sehrwohl Fortran-Compiler, die COMPLEX*16 verarbeiten können, aber Prospero-Fortran versteht in den Versionen bis mmg2.14 nur REAL*8 . Selbst die Zusatz-Library zur Ausnutzung des Coprozessors 68881, die von Prospero angeboten wird, beherrscht nur COMPLEX*8, also COMPLEX in normaler Genauigkeit. Lediglich die Fa. Plünnecke in 3325 Lengede bietet ein Paket von SUBROUTINEn an, das ähnlich meiner Bibliothek komplexe Berechnungen aus- führt. Wer wie ich darauf verzichten will, kann meine Unterprogramme benutzen. Ich liefere Ihnen hier nur den Quellcode (383 Programmzeilen); Sie können ihn in den Editor Ihres Fortran-Systems laden (egal ob ProFortran von Prospero, AC-Fortran von ABSoft oder Swifte-Fortran von GDAT), von dort aus compilen und an Ihr Hauptprogramm linken. Haben Sie Geduld, während der Compiler am Hilfsprogramm CDKETTE hart zu beißen hat! Da die Funktion F vom selben Typ wie das Argument X sein soll, nämlich COMPLEX*16, ist es nicht möglich, FUNCTION-Unterprogramme zu schreiben: COMPLEX*16 FUNCTION F(X) akzeptiert der ProFortran-Compiler nicht. Deshalb sind die Unterprogramme alle SUBROUTINEn. Die doppelt-genaue komplexe Zahl X hat den REAL*8-Realteil A und den REAL*8-Imaginärteil B ; das doppelt-genaue komplexe Ergebnis der Funktionsberechnung F(X) hat den REAL*8-Realteil C und den REAL*8- Imaginärteil D . Diese Variablen müssen als solche in ihrem Typ dekla- riert sein; sie dürfen natürlich auch anders heißen. Dann kann die Funktion F(X) aufgerufen werden mit CALL F(A,B,C,D) Die Eingabewerte A und B werden von der SUBROUTINE nur gelesen und demnach nicht verändert. Doch bedenken Sie, daß REAL-Operationen nie exakt sind, auch wenn sie mit DOUBLE PRECISION ausgeführt werden. Der größte Vorteil von REAL*8 ist der zulässige Wertebereich von ± ( 2.3D-308 ÷ 1.79D308 ) und die exakte 0 . Interpretieren Sie Ergebnisse > ± 2.D305 als ∞ bzw. als vom Rechner nicht mehr korrekt darstellbar. Häufig wird von der komplexen Zahl der Betrag gebildet SQRT(A^2 + B^2) ; dabei darf der Radikand die oben genannten Grenzen des Wertebereichs nicht überschreiten, so daß 1.5D-154 < |A|,|B| ≤ 9.4D153 erfüllt sein muß; s. dazu CDA..., CDDIV, CDLOG, CDSQRT. Ich habe zwar einige Sicherheitsabfragen eingebaut, aber im Interesse einer endlichen Rechenzeit wird nicht alles abgeprüft. Ebenso läßt die Genauigkeit der trigonometrischen INTRINSIC-Funktio- nen zu wünschen übrig, sofern das Argument groß ist. Wenn cos(π/2) = 0 ist, so gilt das für cos(10000π/2) noch lange nicht; eine Glättungs- routine wie CD0D0 hat hier auch keinen Sinn mehr, das Argument vorher modulo 2π zu berechnen erhöht die Genauigkeit ebenfalls nicht. Seien Sie also vorsichtig mit Ergebnissen, bei deren Berechnung trigo- nometrische Funktionen großer Argumente vorkommen! Auch hier ließe sich einiges verbessern, insbesondere durch Reihenentwicklungen am Rand des Wertebereichs. Verwendete Hilfsmittel: Eigenes Gehirnschmalz und der Klassiker Abramowitz, Milton; Irene Stegun u.a. Handbook of mathematical functions Dover Publications, New York, 9.Aufl. 1970 Meine Empfehlung: Kaufen oder benutzen Sie nie die ersten Auflagen solcher Nachschlagewerke; es kostet Stunden, Tage, Wochen, bis Sie merken, daß Sie einem Druckfehler aufgesessen sind! Die von mir geschriebenen Unterprogramme haben die Namen, die diese Funktionen als spezifische Namen üblicherweise haben. Sie als EXTERNAL zu deklarieren ist nicht nötig, weil ProFortran sie ja nicht als INTRINSIC-Funktionen kennt. Die Funktionen im einzelnen, sofern sie komplexe Argumente haben: ================================================================= generischer Name spezifischer Name mein Programmname _______________________________________________________________________ ABS CDABS - Berechnet den doppelt-genauen Betrag Y der COMPLEX*16-Zahl X . Y muß dazu als REAL*8 deklariert worden sein. Ich habe dazu kein Unterprogramm geschaffen, weil sich diese Aufgabe in 1 Zeile erledigen läßt: Y = SQRT(A*A+B*B) - AIMAG - Gibt den Imaginärteil von X als REAL*4-Zahl an; demnach ist AIMAG(X) gleich REAL(B) . CMPLX - - Wandelt die doppelt-genaue komplexe Zahl X in eine COMPLEX*8-Zahl Y um, die natürlich als solche typdeklariert worden sein muß. Auch diese Aufgabe läßt sich in 1 Zeile lösen: Y = CMPLX(A,B) CONJG DCONJG - Verwandelt die komplexe Zahl X = A + i B in die konjugiert- komplexe A - i B . Da bei uns X als Zahlenpaar A,B dargestellt werden muß, ist das Ergebnis der Konjugation A,-B . DBLE - - Wandelt den Realteil der komplexen Größe in eine REAL*8-Zahl; der Imaginärteil wird gleich 0 gesetzt. Wir haben es aber schon mit doppelt-genauen Werten zu tun: also ist DBLE(X) gleich A . DCMPLX - - Diese Funktion habe die zwei Argumente Y und Z , die INTEGER*? oder REAL*? sein können. Sie werden zu REAL*8-Zahlen gemacht, die Realteil und Imaginärteil der entstehenden COMPLEX*16-Größe X dar- stellen. Da wir diese komplexe Variable X aus A und B zusammen- setzen müssen, ist A = DBLE(Y) und B = DBLE(Z) . - DIMAG - Gibt den Imaginärteil von X als REAL*8-Zahl an; demnach ist DIMAG(X) gleich B . - DREAL - Gibt den REAL*8-Realteil der komplexen Zahl an; da sich unsere COMPLEX*16-Zahl X aus den zwei REAL*8-Größen A und B zusammensetzt, ist DREAL(X) gleich A . INT - - Wandelt den Realteil der komplexen Größe in eine INTEGER*4-Zahl. Da der Realteil A unserer Größe X bekannt ist, lautet das Ergebnis: INT(X) gleich INT(A) . REAL REAL - Gibt den REAL*4-Realteil der komplexen Zahl an; da sich unsere COMPLEX*16-Zahl X aus den zwei REAL*8-Größen A und B zusammensetzt, ist REAL(X) gleich REAL(A) . ACOS - CDACOS Standardmäßig gibt es keine Funktion, die den Arcus-Cosinus einer komplexen Größe X = A + i B berechnet, unabhängig davon ob X COMPLEX*8 oder *16 ist. Bei der Programmierung erkannte ich auch warum: Die zugrundegelegte Formel ist numerisch äußerst diffizil. Da SQRT((A+1)^2+B^2) - SQRT((A-1)^2+B^2) durch Rechenfehler nicht gleich 0 werden darf, beschränkt sich der zulässige Wertebereich auf 2.D-8 ≤ |A|,|B| ≤ 1.D15 mit 2.D-8 ≤ |A/B| ≤ 5.D7 . Sofern A aut B gleich 0 ist, lautet die Einschränkung lediglich: 1.D-15 ≤ Betrag des anderen ≤ 1.3D154 (aut ist das ausschlie- ßende Oder). Darüberhinaus darf natürlich X = 0 sein. Jenseits dieser Grenzen ist entweder das Ergebnis falsch oder die Fehler unannehmbar groß. Erwarten Sie deshalb von CDACOS nicht das Unmögliche! Der Realteil C des Ergebnisses ist der Hauptwert des Arcus-Cosinus, also 0 ≤ C ≤ π . Aufruf: CALL CDACOS(A,B,C,D) - - CDACOSH Wie auch CDACOS ist CDACOSH bei Fortran nicht vorgesehen. Zur Berechnung des Area-Cosinus hyperbolicus gilt dasselbe, was ich zu CDACOS geschrieben habe, weil diese Funktion aufgerufen wird. Der Imaginärteil D des Ergebnisses ist der Hauptwert des Arcus- Cosinus, also 0 ≤ D ≤ π . Aufruf: CALL CDACOSH(A,B,C,D) - - CDACOT Die Funktion CDACOT ist bei Fortran gar nicht vorhanden. Zur Berechnung des Arcus-Cotangens gilt dasselbe, was ich zu CDATAN geschrieben habe, weil diese Funktion aufgerufen wird. Der Realteil C des Ergebnisses ist der Hauptwert des Arcus- Cotangens, also - π/2 < C ≤ π/2 . Aufruf: CALL CDACOT(A,B,C,D) - - CDACOTH Wie auch CDACOT ist CDACOTH bei Fortran nicht vorgesehen. Zur Berechnung des Area-Cotangens hyperbolicus gilt dasselbe, was ich zu CDACOT geschrieben habe, weil diese Funktion aufgerufen wird. Der Imaginärteil D des Ergebnisses ist der Hauptwert des Arcus- Cotangens, also - π/2 < D ≤ π/2 . Aufruf: CALL CDACOTH(A,B,C,D) ASIN - CDASIN Wie auch CDACOS ist CDASIN bei Fortran nicht vorgesehen. Zur Berechnung des Arcus-Sinus gilt dasselbe, was ich schon zu CDACOS geschrieben habe, weil beide Funktionen auf dieselbe Berechnungs- routine zurückgreifen. Der Realteil C des Ergebnisses ist der Hauptwert des Arcus-Sinus, also - π/2 ≤ C ≤ π/2 . Aufruf: CALL CDASIN(A,B,C,D) - - CDASINH Wie auch CDASIN ist CDASINH bei Fortran nicht vorgesehen. Zur Berechnung des Area-Sinus hyperbolicus gilt dasselbe, was ich zu CDASIN geschrieben habe, weil diese Funktion aufgerufen wird. Der Imaginärteil D des Ergebnisses ist der Hauptwert des Arcus- Sinus, also - π/2 ≤ D ≤ π/2 . Aufruf: CALL CDASINH(A,B,C,D) ATAN - CDATAN Auch der Arcus-Tangens eines komplexen Arguments ist keine Stan- dardfunktion bei Fortran. In den beiden Polen A = 0 , B = ± 1 liefert CDATAN ± i 1.D308 ; in ihrer Nähe |A| < 1.5D-154 , B = ± 1 entsteht Underflow im Argument des Logarithmus, was dann als Meldung auf dem Bildschirm erscheint. Eine Anregung: Hier könnten Sie durch Umformungen des Logarithmus Abhilfe schaffen. Um Overflow zu umgehen, darf das Quadrat des Argument-Betrags ( |X|^2 ) 1.D308 nicht übertreffen. Der Realteil C des Ergebnisses ist der Hauptwert des Arcus-Tangens, also - π/2 ≤ C ≤ π/2 . Aufruf: CALL CDATAN(A,B,C,D) - - CDATANH Auch CDATANH ist bei Fortran nicht vorgesehen. Zur Berechnung des Area-Tangens hyperbolicus gilt dasselbe, was ich zu CDATAN geschrieben habe, weil diese Funktion aufgerufen wird. Der Imaginärteil D des Ergebnisses ist der Hauptwert des Arcus- Tangens, also - π/2 ≤ D ≤ π/2 . Aufruf: CALL CDATANH(A,B,C,D) COS CDCOS CDCOS Bei der Berechnung des Cosinus komplexer Argumente, deren Realteil in Bogenmaß gegeben sein muß, treten e-Funktionen auf. Deshalb muß |B| ≤ 704 sein; andernfalls tritt zwar kein Overflow ein, ist das Ergebnis aber nicht mehr ganz korrekt (Mehr als e^704 kann der Rechner nun mal nicht richtig). Aufruf: CALL CDCOS(A,B,C,D) COSH - CDCOSH Die Funktion Cosinus hyperbolicus eines komplexen Arguments exi- stiert bei Fortran nicht - hier haben Sie sie! Bedenken Sie, daß der cosh exponentiell wächst und damit sehr schnell sehr groß wird; für |A| trifft hier das zu, was ich unter CDCOS zu |B| geschrieben habe, weil CDCOS zur Berechnung herangezogen wird. Aufruf: CALL CDCOSH(A,B,C,D) - - CDCOT Auch die Cotangens-Funktion kommt bei Fortran nicht vor. Meine Ausführung benutzt das Unterprogramm CDDIV nach Berechnung von cos und sin und liefert die Pole als ± 1.D308 . Das Vorzeichen in diesen ja ungeraden Polen ist dasjenige, das der Cosinus hat. Näheres zu diesem Wert 1.D308 habe ich unter CDTAN ausgeführt. Die hyperbolischen Funktionen ergeben keine Schwierigkeiten, da nur ihr Quotient gebraucht wird, der für große B gegen 1 geht. Aufruf: CALL CDCOT(A,B,C,D) - - CDCOTH Zur Berechnung des Cotangens hyperbolicus eines komplexen Argu- ments wird CDCOT aufgerufen; Details finden Sie dort, wobei B gegen A zu tauschen ist. Aufruf: CALL CDCOTH(A,B,C,D) - - CDDIV Es berechnet den Quotienten der doppelt-genauen komplexen Zahlen A1 + i B1 ( = Zähler) und A2 + i B2 ( = Nenner). Wenn der Nenner gleich 0 ist oder Underflow bei der Betragsbildung auftritt, erscheint eine Meldung, und das Programm gibt den mathematisch falschen Wert C = 0. und D = 0. zurück, sofern Sie "Continue?" mit y beantworten. Und um Overflow zu vermeiden, darf das Quadrat des Nennerbetrags ( A2^2 + B2^2 ) 1.D308 nicht überschreiten. Aufruf: CALL CDDIV(A1,B1,A2,B2,C,D) EXP CDEXP CDEXP Ist die Exponentialfunktion e^X . Sie ist für Argumente endlichen Betrags stets ungleich 0 , und damit das auch numerisch zutrifft, muß |A| ≤ 704 sein. Übertritt A diese Grenzen, wird A := ± 704 gesetzt, damit die Rechnung sicher fortgeführt werden kann; das Ergebnis ist dann aber nicht mehr 100 % richtig. Aufruf: CALL CDEXP(A,B,C,D) LOG CDLOG CDLOG Die Funktion LOG ist der natürliche Logarithmus ln X , also der zur Basis e . Im Komplexen ist der Logarithmus unendlich viel- deutig, so daß mein Programm den Hauptwert ermittelt, d.h. der Imaginärteil D liegt im Bereich - π < D ≤ π . Außerdem ist der Logarithmus für 0 nicht definiert. In dem Fall, oder wenn bei der Betragsbildung Underflow eintritt, unterbricht das Programm die Berechnung, und es erscheint eine Meldung auf dem Bildschirm. Wenn Sie "Continue?" mit y beantworten, liefert CDLOG das mathematisch nicht korrekte Ergebnis C = - 1.D308 und D = 0. (unkorrekt deshalb, weil C = - ∞ sein müßte). Zum Sinn des Werts 1.D308 s. CDTAN. Sofern A ungleich 0 ist (Das Ungleich-Zeichen gibt's in ASCII nicht.), muß das Betragsquadrat des Arguments ( |X|^2 ) ≤ 1.D308 sein, damit kein Overflow eintritt. Aufruf: CALL CDLOG(A,B,C,D) - - CDMUL Dieses Unterprogramm multipliziert die beiden doppelt-genauen komplexen Zahlen A1 + i B1 und A2 + i B2 miteinander. Das Produkt heißt C + i D . Aufruf: CALL CDMUL(A1,B1,A2,B2,C,D) - - CDPOT In doppelter Genauigkeit wird die komplexe Zahl A1 + i B1 ( = Ba- sis) zur komplexen Potenz A2 + i B2 ( = Exponent) erhoben. Auch hier ist Vorsicht geboten, da die Potenzierung sehr schnell den Wertebereich von REAL*8 überschreiten kann. Wenn Basis und Exponent gleich 0 sind, gibt das Programm das mathematisch falsche Ergebnis C = 0. und D = 0. zurück, sofern Sie die Meldung auf dem Bild- schirm nach "Continue?" mit y beantworten. CDPOT berechnet nur 1 Lösung, obwohl eine beliebige Potenz im Kom- plexen beliebig vieldeutig ist: z.B. hat die n. Wurzel n Lösungen. Das ist das Stichwort - mit A2 = 1.D0/n und B2 = 0. kann aus jeder Zahl A1 + i B1 (auch mit B1 = 0. ) die n. Wurzel gezogen werden. Im übrigen gilt, was ich zu CDLOG geschrieben habe, da CDPOT darauf zurückgreift. Aufruf: CALL CDPOT(A1,B1,A2,B2,C,D) SIN CDSIN CDSIN Bei der Berechnung des Sinus eines komplexen Arguments, dessen Realteil in Bogenmaß gegeben sein muß, treten e-Funktionen auf. Deshalb gilt dasselbe, was ich schon zu CDCOS gesagt habe. Aufruf: CALL CDSIN(A,B,C,D) SINH - CDSINH Auch diese Funktion ist bei Fortran nicht für komplexe Argumente eingerichtet. Der Sinus hyperbolicus wächst exponentiell, was den Bereich der Argumente mit richtigen Ergebnissen einengt, s. CDCOSH. Aufruf: CALL CDSINH(A,B,C,D) SQRT CDSQRT CDSQRT Berechnet die doppelt-genaue komplexe Wurzel der COMPLEX*16-Zahl X = A + i B ; die Wurzel heißt C + i D . Da die Quadratwurzel zweideutig ist ( Wurzel(1) = ± 1 ), liefert CDSQRT den Hauptwert, d.h. denjenigen mit C ≥ 0 <=> - π/2 < arctan(D/C) ≤ π/2 . Wenn bei der Betragsbildung Underflow auftritt, erscheint eine Meldung; wenn Sie dann "Continue?" mit y beantworten, wird als Ergebnis C = 0. und D = 0. zurückgegeben. Sofern A ungleich 0 ist, muß |X|^2 ≤ 1.D308 sein, damit kein Overflow auftritt. Aufruf: CALL CDSQRT(A,B,C,D) TAN - CDTAN Die Tangens-Funktion komplexen Arguments ist bei ProFortran nicht vorhanden. Da der Tangens bei den reellen Argumenten A = ± π/2 + k * π , k eine ganze Zahl, und B = 0. die Funktions- werte ± ∞ annimmt und diese Größe ja nicht darstellbar ist, liefert CDTAN die ungefähr größte in DOUBLE PRECISION vorhandene Zahl, nämlich ± 1.D308 . Diese Pole des Tangens sind ungerade, einen exakten Wert haben sie also nicht; entsprechend der Konven- tion gibt CDTAN jedoch 1.D308 mit dem Vorzeichen zurück, das der Sinus gleichen Arguments hat. Die zwei Zahlen ± 1.D308 entsprechen tatsächlich ± ∞ , denn es gilt 1/( ± 1.D308) = ± 1.D-308 ≡ 0. Zur Nicht-Problematik eines großen Imaginärteils B s. CDCOT. Aufruf: CALL CDTAN(A,B,C,D) TANH - CDTANH Mein Tangens hyperbolicus eines komplexen Arguments X ruft CDTAN auf. Es gilt also allgemein, was unter CDTAN, und speziell für große A, was unter CDCOTH zu lesen ist. Aufruf: CALL CDTANH(A,B,C,D) * * *